Вывод уравнений движения ЛА типа "самолет"
Допущения:
1) Конструкция ЛА абсолютно жесткая (ЛА - твердое тело);
2) Атмосфера спокойная (W=0 - скорость ветра);
3) Масса ЛА постоянная т.к. заданное время переходных процессов меньше 60 с;
4) Земля плоская и не вращается;
5) Угол тангажа не равен +/-90 градусов;
6) Спарведлива гипотеза стационарности (согласно этой гипотезе, аэродинамическая сила и аэродиномический момент в нестационарном режиме полета полностью определяются величинами параметров в данный момент времени).
Системы координат и матрицы перехода

MatrixForm[M1={{Cos[Bteta[t]],-Sin[Bteta[t]],0},{Sin[Bteta[t]],Cos[Bteta[t]],0},{0,0,1}}] (*Матрица поворота на угол наклона траектории*)

MatrixForm[M2={{Cos[BPsi[t]],0,Sin[BPsi[t]]},{0,1,0},{-Sin[BPsi[t]],0,Cos[BPsi[t]]}}] (*Матрица поворота на угол пути*)

MatrixForm[M=M2.M1] (*Матрица перехода от траекторной системы координат к нормальной системе координат *)

MatrixForm[ReplaceAll[M,{BPsi[t]->0,Bteta[t]->0}]] (*Проверка подстановкой нулевых углов поворотов. Должна получиться единичная матрица*)

MatrixForm[MM1={{Cos[LPsi[t]],0,Sin[LPsi[t]]},{0,1,0},{-Sin[LPsi[t]],0,Cos[LPsi[t]]}}] (*Матрица поворота на угол рыскания*)

MatrixForm[MM2={{Cos[LTeta[t]],-Sin[LTeta[t]],0},{Sin[LTeta[t]],Cos[LTeta[t]],0},{0,0,1}}](*Матрица поворота на угол тангажа*)

MatrixForm[MM3={{1,0,0},{0,Cos[LGamma[t]],-Sin[LGamma[t]]},{0,Sin[LGamma[t]],Cos[LGamma[t]]}}] (*Матрица поворота на угол крена*)

MatrixForm[MM=MM1.MM2.MM3] (*Матрица перехода от связанной системы координат к нормальной системе координат*)

MatrixForm[SM1={{Cos[Alpha[t]],-Sin[Alpha[t]],0},{Sin[Alpha[t]],Cos[Alpha[t]],0},{0,0,1}}] (*Матрица поворота на угол атаки*)

MatrixForm[SM2={{Cos[Betta[t]],0,Sin[Betta[t]]},{0,1,0},{-Sin[Betta[t]],0,Cos[Betta[t]]}}] (*Матрица поворота на угол скольжения*)

MatrixForm[SM=Transpose[SM2.SM1]] (*Матрица перехода от скоростной системы координат к связанной системе координат*)

Кинематические уравнения поступательного движения

R={L[t],H[t],dZ[t]}(*Радиус вектор центра масс ЛА, L[t] - дальность полета, H[t] - высота полета, dZ[t] - боковое отклонение*)

vg={vgx[t],vdy[t],vgz[t]} (*Вектор земной сокрости центра масс ЛА в нормальной системе координат*)

MatrixForm[D[R,t]]==MatrixForm[vg] (*Кинематическое уравнение движения центра масс самолета в векторной форме*)

MatrixForm[vzem={vk[t],0,0}] (*Вектор земной скорости в траекторной системе координат*)

MatrixForm[vg=M.vzem] (*Проецируем вектор земной скорости на нормальные оси координат*)

MatrixForm[D[R,t]]==MatrixForm[vg] (*Кинематическое уравнение движения центра масс самолета в векторной форме*)

ukpd1=D[R[[1]],t]==vg[[1]]

ukpd2=D[R[[2]],t]==vg[[2]]

ukpd3=D[R[[3]],t]==vg[[3]]

Динамические уравнения поступательного движения

v={vx[t],vy[t],vz[t]} (*Вектор воздушной скорости ЛА в связанной системе координат*)

f={fx,fy,fz} (*Вектор результирующей сил, действующих на ЛА в полете, в связанной системе координат*)

w={wx[t],wy[t],wz[t]} (*Вектор угловой скорости ЛА в связанной системе координат*)

MatrixForm[massa*(D[v,t]+Cross[w,v])]==MatrixForm[f] (*Динамическое уравнение движения центра масс самолета в векторной форме*)

udpd1=(massa*(D[v,t]+Cross[w,v]))[[1]]==f[[1]]

udpd2=(massa*(D[v,t]+Cross[w,v]))[[2]]==f[[2]]

udpd3=(massa*(D[v,t]+Cross[w,v]))[[3]]==f[[3]]

Кинематические уравнения вращательного движения

proekc={D[LGamma[t],t],0,0}+Transpose[MM2.MM3].{0,D[LPsi[t],t],0}+Transpose[MM3].{0,0,D[LTeta[t],t]} (*Проецируем угловые скорости рыскания, тангажа и крена на оси связанной системы координат*)

ukvd1=w[[1]]==proekc[[1]]

ukvd2=w[[2]]==proekc[[2]]

ukvd3=w[[3]]==proekc[[3]]

Динамические уравнения вращательного движения

MatrixForm[TI={{Ixx,-Ixy,-Ixz},{-Iyx,Iyy,-Iyz},{-Izx,-Izy,Izz}}] (*Тензор инерции*)

Ixz=Izx=Iyz=Izy=0 (*Ось OZ перпендикулярна плоскости симметрии XOY*)

Iyx=Ixy

MatrixForm[TI]

MatrixForm[K=TI.w] (*Главный момент количеств движения в связанной системе координат*)

M={Mx,My,Mz} (*Результирующий момент аэродинамических сил и тяги в связанной системе координат(Момент силы тяжести, центробежных и кориолисовых сил инерции вокруг центра масс ЛА равен нулю)*)

lchudvd=D[K,t]+Cross[w,K] (*Полная производная по времени главного вектора момента количеств движения*)

Полная производная главного вектора момента количеств движения равна главному моменту внежних сил, действующих на тело.

udvd1=lchudvd[[1]]==M[[1]]

udvd2=lchudvd[[2]]==M[[2]]

udvd3=lchudvd[[3]]==M[[3]]

MatrixForm[polnsyst={ukpd1,ukpd2,ukpd3,udpd1,udpd2,udpd3,ukvd1,ukvd2,ukvd3,udvd1,udvd2,udvd3}] (*Полная система уравнений движения ЛА*)

Уточним зависимости для сил и моментов

Cy[Alpha[t],Mah[t],ugolstrel,ugolVobr,dzak[t],drv[t]]=D[Cy[Alpha[t],Mah[t],ugolstrel,ugolVobr,dzak[t],drv[t]],Alpha[t]]*(Alpha[t]-Alpha0) (*Коэффициент подъёмной силы*)

lyambda=(l^2)/S (*Относительное удлинение крыла*)

Cx=Cx0+(Cy[Alpha[t],Mah[t],ugolstrel,ugolVobr,dzak[t],drv[t]]^2)/(lyambda*Pi) (*Коэффициент силы лобового сопротивления*)

q=(ro[H[t]]*(vk[t])^2)/2 (*Скоростной напор*)

Maer={kmx[Betta[t],ugolstrel,deler[t],drn[t],wx[t],wy[t],Mah[t]]*q*s*l,kmy[Betta[t],ugolstrel,deler[t],drn[t],wx[t],wy[t],Mah[t]]*q*s*l,kmz[Cy[Alpha[t],Mah[t],ugolstrel,ugolVobr,dzak[t],drv[t]],ugolVobr,drv[t],dzak[t],wz[t],D[Alpha[t],t],Mah[t]]*q*s*ba} (*Аэродинамический момент*)

faer={-Cx*q*s,Cy[Alpha[t],Mah[t],ugolstrel,ugolVobr,dzak[t],drv[t]]*q*s,Cz[Betta[t],drn[t]]*q*s} (*Аэродинамическая сила в скоростной системе координат*)

MatrixForm[faersv=SM.faer](*Аэродинамическая сила в связанной системе координат*)

MatrixForm[Gsv=Transpose[MM].{0,-G,0}] (*Сила тяжести в связанной системе координат*)

Mx=Maer[[1]]

My=Maer[[2]]

Mz=Maer[[3]]-P[dsg[t],vk[t],H[t]]*yp*Cos[fip]

fx=faersv[[1]]+Gsv[[1]]+P[dsg[t],vk[t],H[t]]*Cos[fip]

fy=faersv[[2]]+Gsv[[2]]+P[dsg[t],vk[t],H[t]]*Sin[fip]

fz=faersv[[3]]+Gsv[[3]]

MatrixForm[polnsystut={ukpd1,ukpd2,ukpd3,udpd1,udpd2,udpd3,ukvd1,ukvd2,ukvd3,udvd1,udvd2,udvd3}] (*Уточненная полная система уравнений движения ЛА*)

vektorsost=Flatten[{R,v,w,{LPsi[t],LTeta[t],LGamma[t]}}] (*Вектор состояния системы*)

MatrixForm[SystFK=Flatten[Simplify[Solve[polnsystut,D[vektorsost,t],Sort->False]]]] (*Разрешаем систему относительно производных*)

MatrixForm[systema=Table[SystFK[[i,1]]==SystFK[[i,2]],{i,1,12}]] (*Система уравнений движения ЛА в форме Коши*)

Представим нелинейную нестационарную модель движения ЛА в матричном виде (Overscript[X, ·])[t]=Ax[t]+Bu[t]

MatrixForm[MatrA=Table[D[systema[[i,2]],vektorsost[[j]]],{i,1,12},{j,1,12}]] (*Матрица состояния "свободного" ЛА - матрица A*)

U={drv[t],deler[t],drn[t],dsg[t]} (*Вектор управления U[t]*)

MatrixForm[MatrB=Table[D[systema[[i,2]],U[[j]]],{i,1,12},{j,1,4}]] (*Матрица передачи управления - B*)